06_TGCN

Try using a Targetted Gene Correlation Network to identify the LFY1 and LFY2 networks

Install…

BiocManager::install(c("AnnotationDBi", "GO.db", "preprocessCore", "impute", "rrvgo", "ComplexHeatmap"))
BiocManager::install(c("sva", "GOSim"))
install.packages("MLmetrics")
remotes::install_local("~/Downloads/GOSim_1.40.0.tgz") # get download link from https://bioconductor.org/packages/3.18/bioc/html/GOSim.html
remotes::install_github('juanbot/CoExpNets') 
# remotes::install_github("aliciagp/TGCN") # don't use this, use mine, see below!

use my modified version of TGCN

# remotes::install_local("~/git/TGCN/", force = TRUE)
# or from web:
remotes::install_github("https://github.com/jnmaloof/TGCN", ref = "dev")

Skipping install of 'TGCN' from a github remote, the SHA1 (6985dd91) has not changed since last install.
  Use `force = TRUE` to force installation

library(TGCN)
library(tidyverse)

── Attaching core tidyverse packages ────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     ── Conflicts ──────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(ggpubr)
library(gridExtra)

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

library(gplots)

Attaching package: ‘gplots’

The following object is masked from ‘package:stats’:

    lowess

library(patchwork)

sample.description <- read.csv("../input/sample.description.csv")

load read counts. These should be log2 cpm or Voom transformed or something similar.

lcpm <- read.csv("../output/log2cpm.csv.gz", row.names = 1, check.names = FALSE) 
rownames(lcpm) <- str_remove(rownames(lcpm), fixed(".v2"))
head(lcpm)

dim(lcpm)

[1] 26704    48

Get the GO info and convert to a list.

GOs <- read_delim("../input/Crichardii_676_v2.1.annotation_info.txt")
GOs <- GOs[match(rownames(lcpm), GOs$transcriptName),] %>% drop_na(transcriptName)
gene2GO.list <- GOs %>% 
  select(transcriptName, GO) %>%
  rowwise() %>%
  mutate(GOlist = str_split(GO, pattern = " ")) %>%
  ungroup() %>%
  mutate(GOlist = set_names(GOlist, transcriptName)) %>%
  pull(GOlist)

Looking at the example notebook on the github site, it looks like genes in rows, samples in columns

My genes of interest.

CrLFY1 <- "Ceric.33G031700"

CrLFY2 <- "Ceric.18G076300"

Subset data to build LFY1 and LFY2 networks. Want to find networks around expression of LFY1 or LFY2, so pull those out and also create matrices that don’t have them.

LFY1.expression <- lcpm[str_detect(rownames(lcpm), CrLFY1),] %>%
  unlist() %>%
  as.vector()

input_for_LFY1 <- lcpm[str_detect(rownames(lcpm), CrLFY1, negate = TRUE),]

LFY2.expression <- lcpm[str_detect(rownames(lcpm), CrLFY2),] %>%
  unlist() %>%
  as.vector()

input_for_LFY2 <- lcpm[str_detect(rownames(lcpm), CrLFY2, negate = TRUE),]

LFY1 network

if(dir.exists("../output/TGCN_LFY1"))   system("rm -r ../output/TGCN_LFY1")

if(!dir.exists("../output/TGCN_LFY1"))  dir.create("../output/TGCN_LFY1")

r.lfy1 <- testAllCutoffs(exprData=input_for_LFY1,
                    target=LFY1.expression,
                    covs=NULL,
                    train.split=0.7,
                    nfolds=5,
                    t=10,
                    path="../output/TGCN_LFY1",
                    targetName="LFY1",
                    tissueName="ALL",
                    seed=3333,
                    cutoffs=10:1,
                    n=100, 
                    m=10, 
                    s=10, 
                    minCor=0.3,
                    maxTol=3,
                    save=T,
                    overwrite=T,
                    approach="enrichment", # the approach selected to complete the seed modules
                    report=F,              # if report=T, an automated report will be created
                    gene2GO=gene2GO.list,
                    cellTypeAnnotation=FALSE)

- Step 1: Create a linear regression model for each gene set of hubs based on their ratio of appearance 
Apply LASSO 10 times for feature selection 

Warning: No enrichment can pe performed - there are no feasible GO terms!Warning: No enrichment can pe performed - there are no feasible GO terms!Warning: No enrichment can pe performed - there are no feasible GO terms!

Number of hubs per ratio of appearance
cutoff10  cutoff9  cutoff8  cutoff7  cutoff6  cutoff5  cutoff4  cutoff3  cutoff2  cutoff1 
       1        1        2        4        5       10       11       15       40      150 
Cutoffs selected are 8 7 6 5 4 3 2 1 
Get final model and the cv error for cutoff 8 
Get final model and the cv error for cutoff 7 
Get final model and the cv error for cutoff 6 
Get final model and the cv error for cutoff 5 
Get final model and the cv error for cutoff 4 
Get final model and the cv error for cutoff 3 
Get final model and the cv error for cutoff 2 
Get final model and the cv error for cutoff 1 

- Step 2: represent train and test error per ratio of appearance model

- Step 3: creating a network for each ratio of appearance where the number of hubs is <=30

 Step 3.1: TGCN creation for ratio of appearance 8 

 Step 3.2: TGCN characterization for cutoff 8 

Warning: No enrichment can pe performed - there are no feasible GO terms!

 Step 3.1: TGCN creation for ratio of appearance 7 

 Step 3.2: TGCN characterization for cutoff 7 

 Step 3.1: TGCN creation for ratio of appearance 6 

 Step 3.2: TGCN characterization for cutoff 6 

 Step 3.1: TGCN creation for ratio of appearance 5 

 Step 3.2: TGCN characterization for cutoff 5 

 Step 3.1: TGCN creation for ratio of appearance 4 

 Step 3.2: TGCN characterization for cutoff 4 

 Step 3.1: TGCN creation for ratio of appearance 3 

 Step 3.2: TGCN characterization for cutoff 3 

r.lfy1$selectRatio$nHubs + r.lfy1$selectRatio$stats

Focus on 5 and 6

p <- lapply(r.lfy1$nets, function(cutoff) cutoff$GOenrich$plotStats) 

ggarrange(p$c5 + theme(text=element_text(size=10)), 
          p$c6 + theme(text=element_text(size=10)), 
          ncol=2, nrow=1, common.legend=T, legend="bottom")

r.lfy1$nets$c5$net$moduleSizeSelectionPlot

r.lfy1$nets$c6$net$moduleSizeSelectionPlot

Correlation with trait, I assume

r.lfy1$nets$c5$net$plotCorr

r.lfy1$nets$c6$net$plotCorr

DT::datatable(r.lfy1$nets$c5$net$modules)

DT::datatable(r.lfy1$nets$c6$net$modules)

knitr::include_graphics("../output/TGCN_LFY1/results/LFY1_ALL_c5_TGCN_crossTabPlot.png")

knitr::include_graphics("../output/TGCN_LFY1/results/LFY1_ALL_c6_TGCN_crossTabPlot.png")

grid.arrange(r.lfy1$nets$c5$GOenrich$plotStats, r.lfy1$nets$c5$GOenrich$plotNterms, nrow=2)

grid.arrange(r.lfy1$nets$c6$GOenrich$plotStats, r.lfy1$nets$c6$GOenrich$plotNterms, nrow=2)

LFY2 network

if(dir.exists("../output/TGCN_LFY2"))   system("rm -r ../output/TGCN_LFY2")

if(!dir.exists("../output/TGCN_LFY2"))  dir.create("../output/TGCN_LFY2")

r.lfy2 <- testAllCutoffs(exprData=input_for_LFY2,
                    target=LFY2.expression,
                    covs=NULL,
                    train.split=0.7,
                    nfolds=5,
                    t=10,
                    path="../output/TGCN_LFY2",
                    targetName="LFY2",
                    tissueName="ALL",
                    seed=3333,
                    cutoffs=10:1,
                    n=100, 
                    m=10, 
                    s=10, 
                    minCor=0.3,
                    maxTol=3,
                    save=T,
                    overwrite=T,
                    approach="enrichment", # the approach selected to complete the seed modules
                    report=F,              # if report=T, an automated report will be created
                    gene2GO=gene2GO.list,
                    cellTypeAnnotation=FALSE)

- Step 1: Create a linear regression model for each gene set of hubs based on their ratio of appearance 
Apply LASSO 10 times for feature selection 
Number of hubs per ratio of appearance
cutoff10  cutoff9  cutoff8  cutoff7  cutoff6  cutoff5  cutoff4  cutoff3  cutoff2  cutoff1 
       0        0        1        1        2        2        5       15       39      152 
Cutoffs selected are 6 5 4 3 2 1 
Get final model and the cv error for cutoff 6 
Get final model and the cv error for cutoff 5 
Get final model and the cv error for cutoff 4 
Get final model and the cv error for cutoff 3 
Get final model and the cv error for cutoff 2 
Get final model and the cv error for cutoff 1 

- Step 2: represent train and test error per ratio of appearance model

- Step 3: creating a network for each ratio of appearance where the number of hubs is <=30

 Step 3.1: TGCN creation for ratio of appearance 6 

 Step 3.2: TGCN characterization for cutoff 6 

 Step 3.1: TGCN creation for ratio of appearance 5 

 Step 3.2: TGCN characterization for cutoff 5 

 Step 3.1: TGCN creation for ratio of appearance 4 

 Step 3.2: TGCN characterization for cutoff 4 

 Step 3.1: TGCN creation for ratio of appearance 3 

 Step 3.2: TGCN characterization for cutoff 3 

Focus on 3

p <- lapply(r.lfy2$nets, function(cutoff) cutoff$GOenrich$plotStats) 

ggarrange(p$c3 + theme(text=element_text(size=10)), 
          p$c4 + theme(text=element_text(size=10)), 
          ncol=2, nrow=1, common.legend=T, legend="bottom")

r.lfy2$nets$c3$net$moduleSizeSelectionPlot

Correlation with trait, I assume

r.lfy2$nets$c3$net$plotCorr

DT::datatable(r.lfy2$nets$c3$net$modules)

knitr::include_graphics("../output/TGCN_LFY2/results/LFY2_ALL_c3_TGCN_crossTabPlot.png")

grid.arrange(r.lfy2$nets$c3$GOenrich$plotStats, r.lfy2$nets$c3$GOenrich$plotNterms, nrow=2)

Heat Maps

LFY1

C5

LFY1.c5 <- read_csv("../output/TGCN_LFY1/results/LFY1_ALL_c5_TGCN.csv")

Rows: 160 Columns: 3── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): hubGene, genes
dbl (1): cor
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

LFY1.c5

All genes

rbind(LFY1=LFY1.expression,lcpm[unique(LFY1.c5$genes),]) %>%
  as.matrix() %>%
  heatmap.2(trace="none", cexRow= 0.6, cexCol=0.7, col="bluered", scale="row")

hubgenes

rbind(LFY1=LFY1.expression,lcpm[unique(LFY1.c5$hubGene),]) %>%
  as.matrix() %>%
  heatmap.2(trace="none", cexRow= 0.6, cexCol=0.7, col="bluered", scale="row")

C6

LFY1.c6 <- read_csv("../output/TGCN_LFY1/results/LFY1_ALL_c6_TGCN.csv")

Rows: 70 Columns: 3── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): hubGene, genes
dbl (1): cor
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

LFY1.c6

All genes

hubgenes

rbind(LFY1=LFY1.expression,lcpm[unique(LFY1.c6$hubGene),]) %>%
  as.matrix() %>%
  heatmap.2(trace="none", cexRow= 0.6, cexCol=0.7, col="bluered", scale="row")

LFY2 C3

LFY2.c3 <- read_csv("../output/TGCN_LFY2/results/LFY2_ALL_c3_TGCN.csv")

Rows: 310 Columns: 3── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): hubGene, genes
dbl (1): cor
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

LFY2.c3

All genes

hubgenes

rbind(LFY2=LFY2.expression,lcpm[unique(LFY2.c3$hubGene),]) %>%
  as.matrix() %>%
  heatmap.2(trace="none", cexRow= 0.6, cexCol=0.7, col="bluered", scale="row")